subroutine CalcLocalSteps()
use DataTime, only: nStepLocal, TimeLocal, TimeStepMin, dtGlobal=>dtOut, &
                    TimeGlobal=>dtOut, ttime, & ! TimeGlobal=>ttime
                    TimeCell, dt2Cell, dtFaceCorr, &
                    Balans0,Balans1, eta, dtCell, CurrentTimeCell, TimeEpsilon, nStepCells, &
                    nNextTimeCell, nNextTimeLocal, dnTimeCell, &
                    nNextTimeLocal, nNextTimeCell, nNextTimeCell1, &
                    TimerCRCC_calcChars
use Data_Mesh, only: nFaces, meshCellCenters, nCells, meshFaceCenters, nInternalFaces, &
                     meshFaceOwner, meshFaceNeighbour
use DataCells, only: CellGlobalStepDone, SubStepCellVars, CharsCellPresent, fluxesCellPresent
use Data_Cfoam, only: Rho,U,E,P, RhoNew,Unew,Enew,Pnew,Tnew, nTimeStepAcceleration, nTimeStepAccelerationStart, &
                      epsilon, epsilonSmoothed, propor, RhoNew,Unew,Enew, rank, Psnew,Usnew,Rhosnew, &
                      Ufon, RoFon,Pfon, &
                      OutputFromOneProc
use DataFaces
use DataPerfectGas, only: PressureAllVars
use DataViscosity, only: CalcViscosityOn
use DataPerfectGas, only: cv, Gam=>Gamma
!use TimeMeasure, only: GetHWTimerId,GetHWTick,StartHWTimer,StopHWTimer, outputHWTimer
implicit none

    real(8) :: CellVars1(5), CellVars2(5)
    
    real(8) :: PresFace1,PresCell1,PresFace2,PresCell2,PresFace3
    real(8) :: Flux1(5),Flux2(5),FluxCorrection1(5),FluxCorrection2(5)
    
    real(8) :: minX, maxX

    integer :: iCell,iFace, i,j
    real(8) :: r(3),c, time, Q,dQdT, P0,U0(3)
    
    
    
    integer :: Timer_asyn,Timer_PredStep
    integer :: Timer_Viscous, Timer_CalcChar,Timer_CalcFluxes
    integer :: Timer_CorrStep,Timer_CalcTimeLocal,Timer_transOut,Timer_screenOut
    real :: Secs, cpu_asyn
	integer :: NCalls
	real(8) :: Ticks
	integer(kind=8) :: nfreq
    character(20) :: TimerName

    
!    call GetHWTimerId('Timer_asin',Timer_asyn)
!    call GetHWTimerId('Timer_PredStep',Timer_PredStep)
!    call GetHWTimerId('Timer_Viscous',Timer_Viscous)
!    call GetHWTimerId('Timer_CalcChar',Timer_CalcChar)
!    call GetHWTimerId('Timer_CalcFluxes',Timer_CalcFluxes)
!    call GetHWTimerId('Timer_CorrStep',Timer_CorrStep)
!    call GetHWTimerId('Timer_CalcTimeLocal',Timer_CalcTimeLocal)
!    call GetHWTimerId('Timer_transOut',Timer_transOut)
!    call GetHWTimerId('Timer_screenOut',Timer_screenOut)
!    
!    call GetHWTimerId('TimerCRCC_calcChars',TimerCRCC_calcChars)

!    call StartHWTimer(Timer_asyn)

    CellGlobalStepDone = .false.
    fluxesCellPresent = .true.

!if(allocated(FluxConserv).eq..false.) &
!allocate(FluxConserv(5,2,nInternalFaces))
!FluxConserv = 0d0

!    call OutputTransducerAsynData()
!    call CalcSubStepVars()

    fluxesFacePresent = .false.
    CharsCellPresent = .false.
!    call CalcBalans(Balans0)

!    minX = minVal(abs(meshCellCenters(1,:)))
!    maxX = maxVal(meshCellCenters(1,:))

    do iCell = 1,nCells
        TimeCell(iCell) = 0d0
        nNextTimeCell(iCell) = int(dtCell(iCell)/TimeStepMin + TimeEpsilon)
        dnTimeCell(iCell) = nNextTimeCell(iCell)
    end do
    nNextTimeLocal = minVal(nNextTimeCell(1:nCells))
    !SK nNextTimeCell1(1:nCells) = nNextTimeCell(1:nCells) - nNextTimeLocal ! ==0 = .false./ !=0 = .true.

    do iCell = 1,nCells
        if(nNextTimeCell(iCell) - nNextTimeLocal.eq.0)then
            nNextTimeCell1(iCell) = .false. ! ==0 = .false./ !=0 = .true.
        else
            nNextTimeCell1(iCell) = .true. ! ==0 = .false./ !=0 = .true.
        end if
    end do



!    call StartHWTimer(Timer_asyn)

    do
        
		if(nStepLocal-nTimeStepAccelerationStart.lt.nTimeStepAcceleration) then ! ntimemax
		    propor = max(0d0,1d0*(nStepLocal-nTimeStepAccelerationStart)/nTimeStepAcceleration) ! ntimemax
		    epsilonSmoothed = 1.0d0-propor+propor*epsilon ! utmp(6,100)
	    end if
	    
!        call StartHWTimer(Timer_Viscous)

	    if(CalcViscosityOn)then
            call CalcViscousFluxes()
        end if
	    
!	    call StopHWTimer(Timer_Viscous)
!	    
!	    
!        call StartHWTimer(Timer_PredStep)
    
	    ! predictor-step
	    call CalcSubStepVars()

! ************************************************************* !
            do iCell = 1,nCells
                Tnew(iCell) = pnew(iCell)/(rhonew(iCell)*cv*(Gam-1.0d0))
            end do
            do iFace = 1,nFaces
                psnew(iFace) = FaceVarsNew(5,1,iFace)
                Usnew(1:3,iFace) = FaceVarsNew(2:4,1,iFace)
                Rhosnew(iFace) = FaceVarsNew(1,1,iFace)
            end do
! ************************************************************* !

!        	call StopHWTimer(Timer_PredStep)

!            call MinMaxWrite()

!            if(OutputFromOneProc)then
!                call OutputFromAllProcs()
!            else
!!                write(32) ttime,0,nCells !,(U(1:3,i),p(i),T(i),i=1,nCells) !,E(i),rho(i)
!!                write(32) (U(1:3,i),p(i),T(i),i=1,nCells) !,E(i),rho(i)
!!                write(outRestarts) ttime,nStepInitial,dt,0,nFaces !,(Usnew(1:3,i),psnew(i),Rhosnew(i),i=1,nFaces) !,Esnew(i),Tsnew(i)
!!                write(outRestarts) (Usnew(1:3,i),psnew(i),Rhosnew(i),i=1,nFaces) !,Esnew(i),Tsnew(i)
!!                write(outRestarts) 0,100 !,(utmp(1:8,i),i=1,100)
!!                write(outRestarts) (utmp(1:8,i),i=1,100)
!                call WriteRes()
!            end if
!	    
        ! Calculate Double Set of Characteristics for All Faces
        
!        call StartHWTimer(Timer_CalcChar)

        call CalcChars()
!        call CalcCharsFluxes()
!	    call StopHWTimer(Timer_CalcChar)
        
!        call StopHWTimer(Timer_PredStep)

!	    call StartHWTimer(Timer_CalcFluxes)
        call CalcFluxes() ! Calculate all Flauxess for all faces
!        call StopHWTimer(Timer_CalcFluxes)

        if(CalcViscosityOn)then
            call CalcViscousFluxes()
        end if

!        r(1:3) = meshFaceCenters(1:3,30)
!        write(73,'(8f20.10)') TimeLocal, TimeCell(30),FluxConserv(1,30), TimeCell(31),FluxConserv(2,30), &
!                              abs(FluxConserv(1,30)+FluxConserv(2,30)), &
!                              dt2Cell(31) * fluxCorrection(1,2,30)
!
!        if(fluxesFacePresent(1,30))then
!            write(150,'(2i5,4f20.10)') 1,nStepLocal,TimeLocal,dtCell(30), dtCell(30)*fluxFaces(1,1,30), dtCell(30)*fluxCorrection(1,1,30)
!        end if
!
!        if(fluxesFacePresent(2,30))then
!            write(150,'(2i5,4f20.10)') 2,nStepLocal,TimeLocal,dtCell(31), dtCell(31)*fluxFaces(1,2,30), dtCell(31)*fluxCorrection(1,2,30)
!        end if

!        call CalcBalans(Balans1)

!        call StartHWTimer(Timer_CorrStep)
        ! corrector-step
        call CalcCorrectorStep() !  ! Calculate all variabless in all cells

!        call StopHWTimer(Timer_CorrStep)
!        call CalcBalans(Balans1)
!        write(991,'(6f25.10)') TimeLocal, Balans1(1)/Balans0(1), Balans1(2)/Balans0(2), &
!                               Balans1(3)/Balans0(3), Balans1(4)/Balans0(4), Balans1(5)/Balans0(5)

!        call CalcSubStepVars()

!        call StartHWTimer(Timer_CalcTimeLocal)
        call CalcTimeLocal()
!        call StopHWTimer(Timer_CalcTimeLocal)
        
        nStepLocal = nStepLocal + 1

!        call StartHWTimer(Timer_transOut)
        call OutputTransducerAsynData()
!        call StopHWTimer(Timer_transOut)

!        call StartHWTimer(Timer_screenOut)

        if (mod(nStepLocal,1000).eq.0.and.rank.eq.0) then ! 00
!        if (mod(nStepLocal,1000).eq.0.and.rank.eq.0) then
            call date_and_time(TIME=Eta)
!	        write(*,'(i8,e20.10,a20,2i12)') nStepLocal,ttime+TimeLocal,eta, int(nStepCells/nCells), nNextTimeLocal
	        write(*,'(i8,5e20.10,a20,2i12)') nStepLocal,ttime+TimeLocal,minVal(timeCell(:)),maxVal(timeCell(:)), &
                                    minVal(dtCell(:)),maxVal(dtCell(:)), eta, int(nStepCells/nCells), nNextTimeLocal
!	        write(*,'(i8,e20.10,a20,i8)') nStepLocal,TimeLocal,maxVal(TimeCell(:)),eta, int(nStepCells/nCells)
!	        write(*,'(i8,2e20.10,a20)') nStepLocal,TimeLocal,maxVal(TimeCell(:)),eta
	    end if

!        call StopHWTimer(Timer_screenOut)

!        call StopHWTimer(Timer_CalcFluxes)

!        if (mod(nStepLocal,1000).eq.0) then
!        
!!            call OutputFromAllProcs()

!        
!            write(1,'(f10.5\)',ADVANCE = 'NO') TimeLocal
!            write(2,'(f10.5\)',ADVANCE = 'NO') TimeLocal
!!            do i = 1,10
!!                j = min(1 + i*100,nCells)
!            do j = 1,nCells
!                if(meshFaceCenters(1,j).gt.0.and.abs(meshFaceCenters(2,j)).lt.0.9)then
!                    write(1,'(2f20.10\)',ADVANCE = 'NO') sqrt(sum(meshCellCenters(1:3,j)*meshCellCenters(1:3,j))), Pnew(j)
!                    write(2,'(2f20.10\)',ADVANCE = 'NO') sqrt(sum(meshCellCenters(1:3,j)*meshCellCenters(1:3,j))), Unew(1,j)
!                end if
!            end do
!            write(1,'(a1/)') ' '
!            write(2,'(a1/)') ' '
!
!            write(3,'(f10.5\)',ADVANCE = 'NO') TimeLocal
!            write(4,'(f10.5\)',ADVANCE = 'NO') TimeLocal
!!            do i = 1,10
!!                j = min(775 + i*1000,nFaces)
!            do j = 1,nFaces
!                if(meshFaceCenters(1,j).gt.0.and.abs(meshFaceCenters(2,j)).lt.0.9)then
!                    write(3,'(2f20.10\)',ADVANCE = 'NO') sqrt(sum(meshFaceCenters(1:3,j)*meshFaceCenters(1:3,j))), FaceVarsNew(5,1,j) ! Psnew(j)
!                    write(4,'(2f20.10\)',ADVANCE = 'NO') sqrt(sum(meshFaceCenters(1:3,j)*meshFaceCenters(1:3,j))), FaceVarsNew(2,1,j) ! Usnew(1,j)
!                end if
!            end do
!            write(3,'(a1/)') ' '
!            write(4,'(a1/)') ' '
!        end if

!!	    if (mod(nStepLocal,10).eq.0.and.rank.eq.0) then
!!	        do iCell = 1,nCells
!!                r(1:3) = -meshCellCenters(1:3,iCell)
!!                if(r(1).gt.0d0.and.r(2).gt.0d0.and.r(3).lt.0d0) then
!!                    write(772+nStepLocal/250,'(9f20.10)') TimeCell(iCell),r(1:3), sqrt(sum(r(1:3)*r(1:3))), Pnew(iCell), Unew(1:3,iCell)
!!                end if
!!            end do
!
!!            do iCell = 268,298
!!                r(1:3) = -meshCellCenters(1:3,iCell)
!!                c = sqrt(1.4*PFon/RoFon)
!!                time = TimeCell(iCell) !+dtCell(iCell)
!!                call CalcQ3D(r,time,Q,dQdt,p0,U0)
!!
!!                write(772+iCell,'(10f20.10)') TimeCell(iCell),sqrt(sum(r(1:3)*r(1:3))), &
!!                                              Pnew(iCell), p0, Unew(1:3,iCell), U0(1:3)
!!            end do
!
!!            do iFace = 1,nFaces
!!                r(1:3) = meshFaceCenters(1:3,iFace)
!!                if(r(1).gt.0d0.and.r(2).gt.0d0.and.r(3).lt.0d0) then
!!                    write(772+nStepLocal/250,'(8f20.10)') r(1:3), sqrt(sum(r(1:3)*r(1:3))), Psnew(iFace), Usnew(1:3,iFace)
!!                end if
!!            end do
!!        end if

!        do iCell = 1,nCells !0
!            r(1:3) = -meshCellCenters(1:3,iCell)
!            call CalcQ3D(r,TimeCell(iCell),Q,dQdt,p0,U0) 
!            write(100+nStepLocal,'(7f20.10)') TimeCell(iCell),dtCell(iCell),r(1), Pnew(iCell), Unew(1,iCell), p0,U0(1)
!        end do

!        write(100,'(i9,202f20.10)') nStepLocal,TimeLocal,TimeCell(1:201)
!        write(101,'(i9,202f20.10)') nStepLocal,TimeLocal,dtCell(1:201)
        

!write(772,'(a5,f20.10,i3)') 'Time=', TimeLocal, nStepLocal
!do iFace = 1,nInternalFaces
!    r(1:3) = meshFaceCenters(1:3,iFace)
!    write(772,'(i4,6f20.10)') iFace,r(1), FluxConserv(1:2,iFace), TimeCell(meshFaceOwner(iFace)), &
!              TimeCell(meshFaceNeighbour(iFace)), abs(FluxConserv(1,iFace)+FluxConserv(2,iFace))
!end do
!        if(TimeLocal+TimeStepMin.ge.TimeGlobal) exit ! end local substeps
!        if(TimeLocal+TimeEpsilon.ge.TimeGlobal) exit ! end local substeps
!	    call StopHWTimer(Timer_Viscous)
        if(TimeLocal+TimeEpsilon.ge.dtGlobal) exit ! end local substeps
        
!        call CalcSubStepVars()
    end do
    
!    call StopHWTimer(Timer_asyn)

    
!    call CalcBalans(Balans1)

!    do iFace = 1,nInternalFaces
!        r(1:3) = meshFaceCenters(1:3,iFace)
!        write(772,'(s5,f20.10)') 'Time=', TimeLocal
!        write(772,'(4f20.10)') r(1), FluxConserv(1:2,iFace), abs(FluxConserv(1,iFace)+FluxConserv(2,iFace))
!    end do
    
!    FluxConserv(1,iFace)-FluxConserv(2,iFace)
!    write(772,'(f20.10)') maxVal(abs(FluxConserv(1,:)+FluxConserv(2,:)))

!write(772,'(a5,f20.10,i10)') 'Time=', TimeLocal, nStepLocal
!do iFace = 1,nInternalFaces
!    r(1:3) = meshFaceCenters(1:3,iFace)
!    write(772,'(i4,18f20.10)') iFace,r(1), FluxConserv(1,1:2,iFace), &
!              FluxConserv(2,1:2,iFace), &
!              FluxConserv(3,1:2,iFace), &
!              FluxConserv(4,1:2,iFace), &
!              FluxConserv(5,1:2,iFace), &
!              TimeCell(meshFaceOwner(iFace)), &
!              TimeCell(meshFaceNeighbour(iFace)), abs(FluxConserv(1,1,iFace)+FluxConserv(1,2,iFace)) &
!              , abs(FluxConserv(2,1,iFace)+FluxConserv(2,2,iFace)) &
!              , abs(FluxConserv(3,1,iFace)+FluxConserv(3,2,iFace)) &
!              , abs(FluxConserv(4,1,iFace)+FluxConserv(4,2,iFace)) &
!              , abs(FluxConserv(5,1,iFace)+FluxConserv(5,2,iFace))
!end do

!    CellGlobalStepDone = .false.
!    fluxesCellPresent = .true.
!    call CalcSubStepVars()

! ******************
!    call StartHWTimer(Timer_transOut)
    call OutputTransducerAsynData()
!    call StopHWTimer(Timer_transOut)
! ******************    
!    call StopHWTimer(Timer_asyn)
    
!    write(*,*) 'Summary:'
!    write(*,*) 'Predictor=', Timer_PredStep
!    write(*,*) 'ViscousFlux=',Timer_Viscous
!    write(*,*) 'CalcChars=',Timer_CalcChar
!    write(*,*) 'CalcFluxes=',Timer_CalcFluxes
!    write(*,*) 'Corrector=',Timer_CorrStep
!    write(*,*) 'CalcLocalTime=',Timer_CalcTimeLocal
!    write(*,*) 'TransOut=',Timer_transOut
!    write(*,*) ' '
!    write(*,*) 'Asinchrony=',Timer_asyn


!    call GetHWTick(Timer_asyn,Secs,Ncalls,freq=nfreq,name=TimerName)
!    write(*,'(3a,i8,a,f10.5)') 'Timer ',TimerName,' called ',NCalls, &
!                              ' times CPU time =', secs ! int(100*secs/cpu_dtime), '%'
!    cpu_asyn  = secs
!    do i = 13,21
!        call GetHWTick(i,Secs,Ncalls,freq=nfreq,NTicks=Ticks,name=TimerName)
!        if(ncalls.eq.0) cycle
!        write(*,'(3a,i8,a,f10.5,i5,a)') 'Timer ',TimerName,' called ',NCalls, &
!                                  ' times, CPU time =',Ticks, int(100*Ticks/cpu_asyn), '%'
!    end do
    
!    call outputHWTimer()
    
    write(*,*) ' '

end subroutine CalcLocalSteps

subroutine CalcTimeLocal()
use Data_cfoam, only: nump, cfl
use DataTime, only: TimeLocal, TimeCell
implicit none
    include 'mpif.h'

    if(nump.eq.1) then
        TimeLocal = minVal(TimeCell)
    end if

end subroutine CalcTimeLocal

subroutine CalcSubStepVars()
use Data_Mesh, only: nCells, meshV, meshCells, meshFaceOwner, nInternalFaces
use DataCells, only: sumCellFluxes, fluxesCellPresent, SubStepCellVars, CharsCellPresent, &
                     CellGlobalStepDone
use DataTime, only: dt2Cell, CurrentTimeCell, nStepCells, TimeCell, TimeGlobal=>ttime
use DataFaces, only: ViscousFluxes, CharsFacePresent, fluxesFacePresent, fluxFaces, &
                     FluxConserv, FaceVarsNew
use Data_Cfoam, only: Unew, U, Rhonew, Rho, Enew, E, Pnew, P
use DataPerfectGas, only: PressureAllVars
use DataSource
implicit none

    integer :: iCell, numFace, iFace,iCellOwner, iSide, i, m
    real(8) :: sumViscousCellFluxes(3)
    real(8) :: DivVStress, dtdV
    
    real(8) :: Pi,ek,a,Q, Rho0

    do iCell = 1,nCells
        if(CellGlobalStepDone(iCell)) cycle ! if(fluxesCellPresent(iCell).eq..false.) cycle
        
        if(fluxesCellPresent(iCell).eqv..false.) cycle

!        call CalcViscousFluxesForCell(iCell)

        sumViscousCellFluxes(1:3) = 0d0
        DivVStress = 0d0
        do numFace = 1,6
            iFace = meshCells(numFace,iCell)
            iCellOwner = meshFaceOwner(iFace)
            if(iCellOwner.eq.iCell)then
                iSide = 1
            else
                iSide = 2
            end if
            sumViscousCellFluxes(1:3) = sumViscousCellFluxes(1:3) + ViscousFluxes(1:3,iSide,iFace)

			do m = 1,3
    			DivVStress = DivVStress + FaceVarsNew(m+1,iSide,iFace) * ViscousFluxes(m,iSide,iFace)
			end do

!if(iFace.le.nInternalFaces)then
!    FluxConserv(1:5,iSide,iFace) = FluxConserv(1:5,iSide,iFace) + dt2Cell(iCell) * fluxFaces(1:5,iSide,iFace)
!end if

        end do

        dtdV = dt2Cell(iCell) / meshV(iCell)

        SubStepCellVars(2:5,iCell) = SubStepCellVars(1,iCell) * SubStepCellVars(2:5,iCell) - &
                                     sumCellFluxes(2:5,iCell) * dtdV ! dt2Cell(iCell) / meshV(iCell)
        SubStepCellVars(2:4,iCell) = SubStepCellVars(2:4,iCell) + &
                                     sumViscousCellFluxes(1:3) * dtdV ! dt2Cell(iCell) / meshV(iCell)
        SubStepCellVars(1,iCell) = SubStepCellVars(1,iCell) - &
                                     sumCellFluxes(1,iCell) * dtdV ! dt2Cell(iCell) / meshV(iCell)
        SubStepCellVars(5,iCell) = SubStepCellVars(5,iCell) + DivVStress * dtdV

!	    if(SourcePresent)then
!	        if(iCell.eq.iCellSource)then
!		        call CalcSource(TimeCell(iCell) + dt2Cell(iCell) +TimeGlobal,Q)
!		        Q = Q * dt2Cell(iCell) / meshV(iCell)
!                SubStepCellVars(1,iCell) = SubStepCellVars(1,iCell) + Q
!            end if
!        end if

        SubStepCellVars(2:5,iCell) = SubStepCellVars(2:5,iCell) / SubStepCellVars(1,iCell)

!        CurrentTimeCell(iCell) = CurrentTimeCell(iCell) + dt2Cell(iCell)
        
        nStepCells = nStepCells + 1

!!		if(iCell.eq.1225)then ! 1760)then ! 1875)then
!		if(iCell.eq.28296)then ! 6957)then ! .or.iCell.eq.6958.or.iCell.eq.7117.or.iCell.eq.7118)then ! 1760)then ! 1875)then
!			Pi = 3.1415926535897932384626433832795
!		    ek = 1d7 * 1.739382343d0 / meshV(iCell)
!			a = 10.0d0 
!			Q = ek * sin(2.0d0 * a * Pi * (TimeCell(iCell) + dt2Cell(iCell)))
!			SubStepCellVars(5,iCell) = SubStepCellVars(5,iCell) + Q*dt2Cell(iCell)/SubStepCellVars(1,iCell)
!!			Enew(i) = Enew(i) + Q*dt2/rhonew(i)
!		end if

!	    if(SourcePresent)then
!	        if(iCell.eq.iCellSource)then
!	        	Rho0 = SubStepCellVars(1,iCell)
!!	        	Q = -Q
!	            do i = 2,5
!	                SubStepCellVars(i,iCell) = SubStepCellVars(i,iCell) * Rho0 / (Rho0 - Q)
!                end do
!
!!		        call CalcSource(TimeCell(iCell) + dt2Cell(iCell),Q)
!!		        Q = Q * dt2Cell(iCell)
!!		        SubStepCellVars(1,iCell) = SubStepCellVars(1,iCell) + Q
!!		        SubStepCellVars(2,iCell) = SubStepCellVars(2,iCell) + Q * Unew(1,iCell)
!!		        SubStepCellVars(3,iCell) = SubStepCellVars(3,iCell) + Q * Unew(2,iCell)
!!		        SubStepCellVars(4,iCell) = SubStepCellVars(4,iCell) + Q * Unew(3,iCell)
!!		        SubStepCellVars(5,iCell) = SubStepCellVars(5,iCell) + Q * Enew(iCell)
!
!!                write(18,'(11f20.10)') TimeCell(iCell) + dt2Cell(iCell)+TimeGlobal,Q, Q*SubStepCellVars(2:5,iCell), SubStepCellVars(1:5,iCell) 
!		    end if
!	    end if

	    if(SourcePresent)then
	        if(iCell.eq.iCellSource)then
		        call CalcSource(TimeCell(iCell) + dt2Cell(iCell) +TimeGlobal,Q)
		        Q = Q * dt2Cell(iCell) / meshV(iCell)
                SubStepCellVars(5,iCell) = SubStepCellVars(5,iCell) + Q / SubStepCellVars(1,iCell)
            end if
        end if
        
        Rho(iCell) = RhoNew(iCell)
        U(1:3,iCell) = Unew(1:3,iCell)
        E(iCell) = Enew(iCell)
        P(iCell) = Pnew(iCell)

        RhoNew(iCell) = SubStepCellVars(1,iCell)
        Unew(1:3,iCell) = SubStepCellVars(2:4,iCell)
        Enew(iCell) = SubStepCellVars(5,iCell)
        Pnew(iCell) = PressureAllVars(SubStepCellVars(1:5,iCell))
        
        fluxesCellPresent(iCell) = .false.
    end do

end subroutine CalcSubStepVars

subroutine CalcBalans(Balans)
use DataTime, only: dt2Cell
use Data_Mesh, only: nCells, nFaces, nInternalFaces, meshV, meshFaceOwner
use Data_Cfoam, only: RhoNew,Unew,Enew
use DataFaces, only: fluxFaces, ViscousFluxes
use DataCells, only: fluxesCellPresent
implicit none
    real(8),intent(out) :: Balans(5)

    integer :: iCell,iFace
    real(8) :: dt
    
    Balans(1:5) = 0d0
    do iCell = 1,nCells
        Balans(1) = Balans(1) + RhoNew(iCell) * meshV(iCell)
        Balans(2:4) = Balans(2:4) + RhoNew(iCell) * UNew(1:3,iCell) * meshV(iCell)
        Balans(5) = Balans(5) + RhoNew(iCell) * Enew(iCell) * meshV(iCell)
    end do

!    do iFace = nInternalFaces+1,nFaces
!        iCell = meshFaceOwner(iFace)
!        if(fluxesCellPresent(iCell).eq..false.) cycle
!        dt = dt2Cell(iCell)
!        Balans(1:5) = Balans(1:5) + fluxFaces(1:5,1,iFace) * dt
!    end do

end subroutine  CalcBalans
